#include "../DWMRICore/Point3I.h"
#include "../DWMRICore/Scalar.h"
#include "../DWMRICore/HOT.h"
#include "../DWMRICore/TwoTensor.h"

#include <omp.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

//char *hot_pathname			= "\\\\nabla.cs.purdue.edu\\visdata\\Ziang\\CoherenceStructure\\Datasets\\Alex\\2X\\HOT6.nrrd";
//char *twotensor_pathname	= "\\\\nabla.cs.purdue.edu\\visdata\\Ziang\\CoherenceStructure\\Datasets\\Alex\\2X\\TwoTensor.nhdr";
//char *dst_pathname			= "\\\\nabla.cs.purdue.edu\\visdata\\Ziang\\CoherenceStructure\\Datasets\\Alex\\2X\\Difference.nrrd";
//
//int main()
//{
//	CHOT *pHOT = new CHOT();
//	pHOT->OpenHOTFile(hot_pathname);
//
//	CTwoTensor *pTwoTensor = new CTwoTensor();
//	pTwoTensor->OpenTwoTensorFile(twotensor_pathname);
//
//	uint w = pHOT->GetWidth();
//	uint h = pHOT->GetHeight();
//	uint d = pHOT->GetDepth();
//	uint n = pHOT->GetN();
//
//	float *result = new float[w*h*d];
//	memset(result, 0, sizeof(float)*w*h*d);
//
//	for (uint i = 0; i < w*h*d; ++i) {
//		float tenA[7], tenB[7];
//		pTwoTensor->GetTensorA(i, tenA);
//		pTwoTensor->GetTensorB(i, tenB);
//		if (tenB[0] > 0.01f) {
//			float *hot = new float[n];
//			pHOT->GetHOTFast(i, hot);
//			CVector3F hotDirs[3];
//			if (pHOT->HOT2Directions(hot, hotDirs) >= 2) {
//				CVector3F tenDirs[2];
//				pTwoTensor->Tensor2Direction(tenA, &(tenDirs[0]));
//				pTwoTensor->Tensor2Direction(tenB, &(tenDirs[1]));
//
//				float dif00 = fabs(InnerProduct(hotDirs[0], tenDirs[0]));
//				float dif01 = fabs(InnerProduct(hotDirs[0], tenDirs[1]));
//
//				float dif10 = fabs(InnerProduct(hotDirs[1], tenDirs[0]));
//				float dif11 = fabs(InnerProduct(hotDirs[1], tenDirs[1]));
//
//				float dif0 = dif00 > dif01 ? dif00 : dif01;
//				float dif1 = dif10 > dif11 ? dif10 : dif11;
//
//				float dif = dif0 > dif1 ? dif1 : dif0;
//
//				result[i] = dif;
//			}
//			SafeDeleteArray(hot);
//		}
//	}
//
//	CScalar *pResult = new CScalar();
//	pResult->CreateScalar(w, h, d, result);
//	pResult->SaveScalarFile(dst_pathname);
//
//	SafeDelete(pHOT);
//	SafeDelete(pTwoTensor);
//	SafeDelete(pResult);
//	SafeDeleteArray(result);
//
//	return 0;
//}

//char *src_pathname = NULL;
//
//void PrintUsage()
//{
//	printf("Usage: TDStatistics -src src_file\n");
//}
//
//int ParseArguments(const int argc, char *argv[])
//{
//	for (int i = 1; i < argc; ++i) {
//		if (strcmp(argv[i], "-src") == 0) {
//			src_pathname = argv[++i];
//		} else {
//			return -1;
//		}
//	}
//
//	if (src_pathname == NULL)
//		return -2;
//
//	return 0;
//}
//
//int main(int argc, char *argv[])
//{
//	if (ParseArguments(argc, argv) != 0) {
//		PrintUsage();
//		return 0;
//	}
//
//	CScalar *pSrc = new CScalar();
//	pSrc->OpenScalarFile(src_pathname);
//
//	uint w = pSrc->GetWidth();
//	uint h = pSrc->GetHeight();
//	uint d = pSrc->GetDepth();
//
//	int count = 0;
//	for (uint i = 0; i < w * h *d; ++i) {
//		if (pSrc->GetScalarFast(i) > 0.5f) {
//			count++;
//		}
//	}
//
//	printf("=========== Statistics ===========\n");
//	printf("Total Voxels = %d\n", w*h*d);
//	printf("Avilable Voxels = %d\n", count);
//	printf("Ratio = %.3f%%\n", (float)count/(float)(w*h*d)*100.0f);
//	printf("=============== End ==============\n");
//
//	SafeDelete(pSrc);
//}


/**********************************************************/
/**********************************************************/

char *hot_pathname = NULL;
char *mask_pathname = NULL;
char *dst_pathname = NULL;

void PrintUsage()
{
	printf("Usage: TDStatistics -hot hot_file -mask mask_file -dst dst_file\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-hot") == 0) {
			hot_pathname = argv[++i];
		} else if (strcmp(argv[i], "-mask") == 0) {
			mask_pathname = argv[++i];
		} else if (strcmp(argv[i], "-dst") == 0) {
			dst_pathname = argv[++i];
		} else {
			return -1;
		}
	}

	if (hot_pathname == NULL || mask_pathname == NULL || dst_pathname == NULL)
		return -2;

	return 0;
}


int main(int argc, char *argv[])
{
	if (ParseArguments(argc, argv) != 0) {
		PrintUsage();
		return 0;
	}

	CHOT *pHOT = new CHOT();
	pHOT->OpenHOTFile(hot_pathname);

	CScalar *pMASK = new CScalar();
	pMASK->OpenScalarFile(mask_pathname);

	uint w = pHOT->GetWidth();
	uint h = pHOT->GetHeight();
	uint d = pHOT->GetDepth();
	uint n = pHOT->GetN();

	CScalar *pResult = new CScalar();
	pResult->CreateScalar(w, h, d);

	printf("\n");
	for (uint z = 0; z < d; ++z) {
		for (uint y = 0; y < h; ++y) {
			omp_set_num_threads(16);
#pragma omp parallel for
			for (int x = 0; x < w; ++x) {
				if (pMASK->GetScalarFast(x, y, z) < 0.5f) {
					pResult->SetScalar(x, y, z, 0.0f);
				} else {
					float *hot = new float[n];
					pHOT->GetHOTFast(x, y, z, hot);
					CVector3F dirs[3];
					int count = pHOT->HOT2Directions(hot, dirs);
					if (count == 3) {
						pResult->SetScalar(x, y, z, 1.0f);
					} else {
						pResult->SetScalar(x, y, z, 0.0f);
					}
					//pResult->SetScalar(x, y, z, (float)count);
					delete[] hot;
				}
			}
		}
		float ff = float(z+1) / (float)d * 100.0f;
		printf("\r%.2f%% done.", ff);
		fflush(stdout);
	}
	printf("\r100.0%% done.\n\n");

	pResult->SaveScalarFile(dst_pathname);

	SafeDelete(pHOT);
	SafeDelete(pMASK);
	SafeDelete(pResult);
}